Mouse Gastrulation scATAC-seq & scRNA-seq

## Load libraries


library(ArchR)
library(knitr)
library(rhdf5)
library(uwot)
library(tidyverse)
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1
library(scater)
library(zellkonverter)
#library(caret)
h5disableFileLocking()

The scATAC-seq dataset contains two replicates each from 4 stages of mouse gastrulation, E7.5, E8.0, E8.5 and E8.75. We will use the mm10 Genome, because that is the genome.

# names of gastrulation stages
stages <- c("E7.5_rep1", "E7.5_rep2", "E8.0_rep1", "E8.0_rep2", 
 "E8.5_rep1", "E8.5_rep2", "E8.75_rep1", "E8.75_rep2")

# combine all raw fragment files
inputFiles <- setNames(c("raw_files/E7.5_rep1_atac_fragments.tsv.gz",
                       "raw_files/E7.5_rep2_atac_fragments.tsv.gz",
                       "raw_files/E8.0_rep1_atac_fragments.tsv.gz",
                       "raw_files/E8.0_rep2_atac_fragments.tsv.gz",
                       "raw_files/E8.5_rep1_atac_fragments.tsv.gz",
                       "raw_files/E8.5_rep2_atac_fragments.tsv.gz",
                       "raw_files/E8.75_rep1_atac_fragments.tsv.gz",
                       "raw_files/E8.75_rep2_atac_fragments.tsv.gz"),
                       c("E7.5_rep1", "E7.5_rep2", "E8.0_rep1", "E8.0_rep2", 
                         "E8.5_rep1", "E8.5_rep2", "E8.75_rep1", "E8.75_rep2"))

addArchRGenome("mm10")

Create Arrow Files

We will read in the accessible fragments from the fragments file. QC information is automatically computed for each cell. Based on visual inspection of TSS Enrichment score vs log10(nFrags) I used a threshold of at least 1000 unique fragments and a TSS Enrichment score of at least 2.75. At the beginning it is better to be less stringent, because we can still filter the data at a later stage. The cells are filtered based on these threshold and data are stored in arrow files. Genome-wide tile matrices are created using 500-bp bins. The tile matrix T is binary matrix and contains information whether in a particular 500-bp bin/tile i there were any insertions in a particular cell j.

\(T \in \mathbb{R}^{NXD}\) with \(T_{i,j} = 1\) if an insertion is present in tile \(i\) and cell \(j\). \(T_{i,j} = 0\) if no insertion is present in that cell, meaning \(T\) is a binary matrix.

In addition a GeneScoreMatrix is created.

I set the size of the maximum number of fragments included to 1e+07, instead of the default 1e+05, because we cut off some cells otherwise.

ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = names(inputFiles),
  minTSS =2.75, #Dont set this too high because you can always increase later
  minFrags = 1000, # minimum number of mapped fragments required
  # count matrix, instead of using peak it uses fixed-width sliding window of bins across the whole genome
  addTileMat = TRUE, 
  addGeneScoreMat = TRUE, # uses signal proximal to the TSS to estimate gene activity
  subThreading = FALSE,
  maxFrags = 1e+07,
  minFragSize = 10,
  maxFragSize = 2000,
  QCDir = "QualityControl",
  force = FALSE, # change to TRUE if you want to ovewrite the files
  # the length in bp that wraps around nucleosomes -> 
    #identify fragments as sub-nucleosome spanning, mono-nucleosome spanning or multi-nucleosome spanning
  nucLength = 147, 
  # integer vector -> define region up/downstream of TSS to include as promoter region
  # can be used to calculate e.g fraction of reads in promoter (FIP)
  promoterRegion = c(2000, 100),
  # parameters for computing TSS enrichment scores, window (bp) centered at TSS = 101
  # flanking window = 2000 bp
  # norm = size of flank window used for normalization = 100 bp
  # accessibility within 101 bp surrounding the TSS will be normalized to accessibility
  # in 100 bp bins from -2000:-1901bp and 1901: 2000
  TSSParams = list(101, 2000, 100),
  # which chromosomes to exclude form downstream analysis
  # in human and mouse: mitochondrial DNA (chrM) and male sex chromosome (chrY)
  # the fragments are still stored in the arrow files
  excludeChr = c("chrM", "chrY"),
  # number of chunks to divide chromosomes in -> low-memory parallelized reading of input files
  nChunk = 5,
  # name of field in input bam file containing the barcode tag information
  bcTag = "qname",
  offsetPlus = 4, # offset applied to + stranded Tn5 insertion -> account for precise Tn5 binding site
  offsetMinus = -5, 
  logFile = createLogFile("createArrows")
)

Inferring Doublets

Should be removed, because they can interfere with downstream analysis.

Doublet detection-and-removal algorithm: Heterotypic doublets are identified by generating a collection of synthetic doublets. These synthetic doublets are projected into low-dimensional embeddings. Searching for nearest neighbours to the synthetic doublets we can identify doublets in the dataset. This outperforms the prediction of doublets using fragment number (ROC-AUC). (Compared to demuxlet as ground truth)

We can also identify doublets in the scRNA-seq space if we have paired data and remove the cells in this way.

# for each sample provided doublet information will be assigned to each cell
# this way we can remove doublet-based clusters downstream
doubScores <- addDoubletScores(
  useMatrix = "TileMatrix",
  input = ArrowFiles,
  k = 10, #Refers to how many cells near a "pseudo-doublet" to count.
  nTrial = 5, # number of time to simulate nCell doublets 
  knnMethod = "UMAP", #Refers to the dimensionality reduciton method to use for nearest neighbor search.
  LSIMethod = 1 # oder of normalization: tf-log(idf)
)

Create ArchRProject

An ArchR Project is initialized with some important attributes:

• ouput directory
• sample names
• sampleColData -> matrix containing data for each sample
• cellColData -> contains data associated with each cell
  • after using addDoubletScore() there will be a column for “Doublet Enrichment” and “Doublet Score”
• total number of cells (excluding doublets)
• median TSS score & median number of fragments across all cells and samples

proj <- ArchRProject(
  ArrowFiles = ArrowFiles, 
  outputDirectory = "gastrArchROutput",
  copyArrows = TRUE, #This is recommened so that you maintain an unaltered copy for later usage.
  geneAnnotation = getGeneAnnotation(),
  #genomeAnnotation = getGeneAnnotation(),
  showLogo = FALSE
)

proj

# read in existing project, so we don't have to rerun the creation
#saveArchRProject(proj, outputDirectory = "01_FirstQC_ArchRobject/", load = FALSE)

proj <- loadArchRProject(path = "01_FirstQC_ArchRobject/")

Quality Control

Have a look at this for additional QC! https://bioconductor.org/packages/devel/bioc/vignettes/ATACseqQC/inst/doc/ATACseqQC.html

1. the number of unique nuclear fragments (as opposed to mitochondrial fragments) A cell with very few usable fragments will not provide enough data to make useful conclusions.
2. signal-to-background ratio: if this is low this probably corresponds to dying cells where the entire genome allows random transposition
3. fragment size distribution: since 147 bp are wrapped around a nucleosome it is expected that there are depletions of fragments of this length at regular intervals. We expect to see a periodic distribution of fragmetn size corresponding to nucleosomes (mono, di, tri, …), because Tn5 cannot cut DNA that is tightly wrapped around a nucleosome.

Data before QC and corresponding plots are saved in the Quality Control output folder.

log10(unique fragments) vs TSS enrichment

before QC

We used a TSS enrichment score of > 2.75 and a number of unique fragments > 1000 (log10(1000) = 3). It might make sense to use different thresholds for different samples.

dir <- "QualityControl/"
#all_files <- list.files(path)


beforeQC <- map(stages, function(n) {
  path <- paste0(dir, n, "/")
  all_files <- list.files(path)
  file <- all_files[grepl(paste0(n, "-Pre-Filter-Metadata.rds"), all_files)]
  # read in file for sample n containing original data before QC
  df <- as.data.frame(readRDS(paste0(dir, n, "/", file)))
  df %>% mutate(Sample = n) %>% 
    ggplot(aes(x = log10(nFrags), y = TSSEnrichment)) +
    geom_density2d_filled() +
    geom_hline(yintercept = 2.75, color="green", linetype="dashed") +
    geom_vline(xintercept = log10(1000), color="green", linetype="dashed") +
    #geom_xsidedensity(aes(x=log10(pre_filter_meta$nFrags))) +
    #geom_ysidedensity(aes(y = pre_filter_meta$TSSEnrichment)) +
    ylim(c(2.5, 4.5)) +
    facet_wrap(~Sample) +
    theme(legend.position = "none") +
    labs(x = "Log10 Unique Fragments", y = "TSS Enrichment Score")
})
  
gridExtra::grid.arrange(grobs = beforeQC, ncol = 2)
## Warning: Removed 361 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 382 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 124 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 103 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 75 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 62 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 269 rows containing non-finite values (stat_density2d_filled).
## Warning: Removed 277 rows containing non-finite values (stat_density2d_filled).

after QC

metadata <- as.data.frame(getCellColData(proj))
metadata %>% head
##                                 Sample TSSEnrichment ReadsInTSS ReadsInPromoter
## E7.5_rep2#TGCTGGATCAAGCTAC-1 E7.5_rep2         3.428      32719          132988
## E7.5_rep2#GCCACTAAGTGAGCAA-1 E7.5_rep2         2.759      12761           61758
## E7.5_rep2#ATTACGTCATGGCCCA-1 E7.5_rep2         2.945      13192           60823
## E7.5_rep2#GATTCAGGTAACTACG-1 E7.5_rep2         3.365      15631           64299
## E7.5_rep2#GTGCTGATCCCTCGCA-1 E7.5_rep2         2.810       6218           30111
## E7.5_rep2#CTGCAATAGTAACGAG-1 E7.5_rep2         2.996       7454           33941
##                              ReadsInBlacklist PromoterRatio PassQC
## E7.5_rep2#TGCTGGATCAAGCTAC-1            17749    0.15421294      1
## E7.5_rep2#GCCACTAAGTGAGCAA-1            14089    0.07720696      1
## E7.5_rep2#ATTACGTCATGGCCCA-1            14215    0.09650094      1
## E7.5_rep2#GATTCAGGTAACTACG-1             7636    0.15642090      1
## E7.5_rep2#GTGCTGATCCCTCGCA-1             8230    0.07454879      1
## E7.5_rep2#CTGCAATAGTAACGAG-1             7397    0.09159430      1
##                              NucleosomeRatio nMultiFrags nMonoFrags nFrags
## E7.5_rep2#TGCTGGATCAAGCTAC-1        1.547911       89468     169230 431183
## E7.5_rep2#GCCACTAAGTGAGCAA-1        2.864000       89587     103507 399951
## E7.5_rep2#ATTACGTCATGGCCCA-1        1.202697       59343     143071 315142
## E7.5_rep2#GATTCAGGTAACTACG-1        1.288291       40144      89819 205532
## E7.5_rep2#GTGCTGATCCCTCGCA-1        1.069593       35346      97582 201955
## E7.5_rep2#CTGCAATAGTAACGAG-1        1.177268       32949      85097 185279
##                              nDiFrags DoubletScore DoubletEnrichment
## E7.5_rep2#TGCTGGATCAAGCTAC-1   172485      0.00000         0.1000000
## E7.5_rep2#GCCACTAAGTGAGCAA-1   206857     19.43591         2.8166667
## E7.5_rep2#ATTACGTCATGGCCCA-1   112728      0.00000         0.7333333
## E7.5_rep2#GATTCAGGTAACTACG-1    75569      0.00000         0.2333333
## E7.5_rep2#GTGCTGATCCCTCGCA-1    69027    132.05092         6.5000000
## E7.5_rep2#CTGCAATAGTAACGAG-1    67233     69.25591         4.7000000
##                              BlacklistRatio
## E7.5_rep2#TGCTGGATCAAGCTAC-1     0.02058175
## E7.5_rep2#GCCACTAAGTGAGCAA-1     0.01761341
## E7.5_rep2#ATTACGTCATGGCCCA-1     0.02255333
## E7.5_rep2#GATTCAGGTAACTACG-1     0.01857618
## E7.5_rep2#GTGCTGATCCCTCGCA-1     0.02037583
## E7.5_rep2#CTGCAATAGTAACGAG-1     0.01996179

metadata %>% group_by(Sample) %>% summarise(n = n())
## # A tibble: 8 × 2
##   Sample         n
##   <chr>      <int>
## 1 E7.5_rep1   8708
## 2 E7.5_rep2  12919
## 3 E8.0_rep1   7024
## 4 E8.0_rep2   5943
## 5 E8.5_rep1  10955
## 6 E8.5_rep2  11091
## 7 E8.75_rep1  4002
## 8 E8.75_rep2  5553

metadata %>% ggplot() +
  geom_density2d_filled(aes(x = log10(nFrags), y = TSSEnrichment)) +
  geom_hline(yintercept = 2.75, color="green", linetype="dashed") +
  geom_vline(xintercept = log10(1000), color="green", linetype="dashed") +
  #geom_xsidedensity(aes(x=log10(pre_filter_meta$nFrags))) +
  #geom_ysidedensity(aes(y = pre_filter_meta$TSSEnrichment)) +
  ylim(c(2.5, 4.5)) +
  facet_wrap(~Sample) +
  theme(legend.position = "none") +
  labs(x = "Log10 Unique Fragments", y = "TSS Enrichment Score")
## Warning: Removed 246 rows containing non-finite values (stat_density2d_filled).

ArchR plots

```#{r, fig.width=8, fig.height=15} #create 3 separate dataframes for all samples samples <- map(unique(proj\(Sample), function(name){ index <- BiocGenerics::which(proj\)Sample %in% name) cells <- proj$cellNames[index] sample_subset <- proj[cells] df <- getCellColData(sample_subset, select = c(“log10(nFrags)”, “TSSEnrichment”)) p <- ggPoint( x = df[, 1], y = df[, 2], colorDensity = TRUE, # should the density of points on the plot be indicated by color? continuousSet = “sambaNight”, xlabel = “Log10 unique fragments”, ylabel = “TSS enrichment”, title = paste0(“Sample:”, name), #xlim = c(log10(500), quantile(df[,1], probs = 0.99)), #ylim = c(0, quantile(df[,2], probs = 0.99)) ) + geom_hline(yintercept = 2.75, lty = “dashed”) + geom_vline(xintercept = 3, lty = “dashed”) })

#do.call(“grid.arrange”, c(samples, ncol=4))

gridExtra::grid.arrange(grobs = samples, ncol = 2)

### Plotting sample statistics

* when we have distinct samples, it can be important to compare various
metric between samples
* ridge plots & violin plots are used for grouped data


```r

p1 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_")) %>% 
  ggplot() +
  geom_density(aes(x = TSSEnrichment, fill = Sample), alpha = 0.8) 
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.

p2 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_")) %>% 
  ggplot() +
  ggridges::geom_density_ridges(aes(x = TSSEnrichment, y = Sample,
                                    fill = Sample), alpha = 0.8) +
  theme(legend.position = "none")

p3 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_")) %>% 
  ggplot() +
  geom_violin(aes(x = Sample, y = TSSEnrichment, fill = Sample), alpha = 0.8) +
  geom_boxplot(aes(x = Sample, y = TSSEnrichment,fill = Sample), alpha = 0.4) + 
  theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") +
  labs(title = "TSS Enrichment")


cowplot::plot_grid(p3, p2, p1, ncol = 3)

p1 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_"), log10_nFrags = log10(nFrags)) %>% 
  ggplot() +
  geom_density(aes(x = log10_nFrags, fill = Sample), alpha = 0.8)

p2 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_"), log10_nFrags = log10(nFrags)) %>% 
  ggplot() +
  ggridges::geom_density_ridges(aes(x = log10_nFrags, y = Sample,
                                    fill = Sample), alpha = 0.8) +
  theme(legend.position = "none") 


p3 <- as_data_frame(getCellColData(proj)) %>% 
  mutate(Sample = str_remove(Sample, "scATAC_"), log10_nFrags = log10(nFrags)) %>% 
  ggplot() +
  geom_violin(aes(x = Sample, y = log10_nFrags, fill = Sample), alpha = 0.8) +
  geom_boxplot(aes(x = Sample, y = log10_nFrags,fill = Sample), alpha = 0.4) + 
  theme(legend.position = "none") +
  labs(title = "number of fragments") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") 

cowplot::plot_grid(p3, p2, p1, ncol = 3)

Plot Fragment Size Distribution & TSS Enrichment Profiles

• the distribution of fragments size can be very different between samples, cell types and batches -> these differences do not necessarily correlate with differences in quality
• the dip is the fragment size of a nucleosome ~147bp
• TSS enrichment profiles
  • clear peak in the center
  • smaller shoulder peak right of the center caused by well positioned +1 nucleosome

#{r, fig.width=5, fig.height=5} p1 <- plotFragmentSizes(ArchRProj = proj) p2 <- plotTSSEnrichment(ArchRProj = proj) ggAlignPlots(p1, p2, type = "v")

Filtering Doublets

In the plots below the Doublet Enrichment represents the enrichment of simulated dobulets nearby each cell, compared to the expected number of doublets if a uniform distriution is assumed.

The doublet density vissualizes where the synthetic doublets are located after projecting them into the 2D embedding.

# the doublet information is saved in a simpleListobject



dir <- "QualityControl/"


doublet_enrichment <- map(stages, function(n) {
  path <- paste0(dir, n, "/")
  all_files <- list.files(path)
  file <- all_files[grepl(paste0(n, "-Doublet-Summary.rds"), all_files)]
  # read in file containing doublet summary for sample n 
  doublet_summary <- readRDS(paste0(dir, n, "/", file))
  doublet_summary[[1]] %>% 
  ggplot() +
  geom_point(aes(x = X1, y = X2, col = enrichment), size = .8) +
  scale_color_viridis_c() +
  guides(fill=guide_legend(title="DoubletEnrichment")) +
  labs(title = paste0("Simulated Doublet Enrichment", n))
})
  
gridExtra::grid.arrange(grobs = doublet_enrichment, ncol = 2)

doublet_density <- map(stages, function(n) {
  path <- paste0(dir, n, "/")
  all_files <- list.files(path)
  file <- all_files[grepl(paste0(n, "-Doublet-Summary.rds"), all_files)]
  # read in file containing doublet summary for sample n 
  doublet_summary <- readRDS(paste0(dir, n, "/", file))
  doublet_summary[[2]] %>% 
  ggplot() +
  geom_point(aes(x = x, y = y, col = density), size = .1) +
  geom_point(data = doublet_summary[[1]], aes(x = X1, y = X2), size = .1,
             alpha = .3) +
  scale_color_viridis_c() +
  guides(fill=guide_legend(title="Simulated Doublet Density")) +
  labs(title = "Simulated doublet density overlayed")

})

gridExtra::grid.arrange(grobs = doublet_density, ncol = 2)

With the function addDoubleScores() information on predicted doublets has been added. Now we will filter the putative doublets. ArchR automatically prints the number of cells removed from each sample and the corresponding percentage.

arguments:

• cutEnrich = minimum cutoff for DoubletEnrichment, number of simulated doublets divided by expected number given a random uniform distribution
• cutScore = minimum cutoff for Doublet Score, represents -log10(binomial adjusted p-value) for the DoubletEnrichmentadd
• filterRatio = maximum ratio of predicted doublets to filter based on number of pass-filter cells (A higher filterRatio means that more cells are removed) e.g. 5000 cells

maximum would be filterRatio * 5000 / 100000 = filterRatio * 5000 * 0.05

This way samples with different percentage of doublets will be filtered accordingly.

As can be seen below, in E7.5_rep1 no doublets are filtered. Why might that be?

# in our case we now have 10 251 cells as opposed to 10 661 cells before
# filtering -> 410 cells were removed (3.85%)
proj <- filterDoublets(ArchRProj = proj)

Use already annotated cells

Subset cells

What is the overlap between my QC and the processed and annotated data?

In total after QC and doublet filtering the dataset contains 59,113 cells. First, we have to remove the cells from E8.75_rep1/2, because these samples are not included in the processed and annotated data. 8,880 cells are from samples E8.75_rep1/2. In the processed data matrix we find 38,158 cells. Only 30,803 cells pass RNA QC. Of the scRNA-seq dataset only 30,282 out of 36,822 cells pass the scATAC-seq QC. This means that we have 30,282 cells with annotated celltypes. These are the cells we want to use in the following analysis. We will, therefore, subset our dataset. We are now left with 28,154 cells for further analysis. Apparently we filtered some cells not fitlered in the processed data.

To visualize this in a better way create a Venn diagram.

# read in processed scATAC-seq data
proc_atac <- readH5AD("anndata_atac.h5ad")

# remove cells which do not pass the rnaQC
remove <- (as.data.frame(colData(proc_atac)) %>% rownames_to_column("cell") %>% 
                          filter(pass_rnaQC == TRUE))[["cell"]]

filt_processed_atac <- proc_atac[ , colnames(proc_atac) %in% remove]

# check the RNA dataset as well to see if the cells overlap with the ATAC cells
proc_rna <- readH5AD("anndata_rna.h5ad")
## Warning in .extract_or_skip_assay(skip_assays = skip_assays, hdf5_backed =
## hdf5_backed, : 'spliced' layer matrix does not support transposition and has
## been skipped
## Warning in .extract_or_skip_assay(skip_assays = skip_assays, hdf5_backed =
## hdf5_backed, : 'unspliced' layer matrix does not support transposition and has
## been skipped


# keep the cells which pass both rnaQC and atacQC
keep <- colnames(proc_rna[,colnames(proc_rna) %in% colnames(proc_atac)])

# remove the samples E8.75
idxSample <- BiocGenerics::which(proj$Sample %in% c("E7.5_rep1", "E7.5_rep2",
                                              "E8.0_rep1", "E8.0_rep2", 
                                              "E8.5_rep1", "E8.5_rep2"))
filt_proj <- proj[proj$cellNames[idxSample]]

In the following plot you can see the cells which passed quality control for scRNA-seq, scATAC-seq and which are included in the ArchR object with a very loose quality control. We will proceed only using the overlap between the three, because for the overlap the celltypes are annotated already. These are also the cells which are used as an input to scDori.

# plot the overlap between the different versions
p <- ggVennDiagram::ggVennDiagram(list(colnames(proc_atac), colnames(proc_rna), filt_proj$cellNames),
              label_alpha = 0,
              category.names= c("processed ATAC", "processed RNA", "ArchR project"),
              #set.size = 20,
              label = "count",
              label_size = 6) +
   scale_fill_distiller(palette = "RdBu") +
  labs(title = "All TFs with different activity")
p

# keep only the cells which passed RNA QC and ATAC QC -> celltype annotations
overlap <- filt_proj$cellNames[filt_proj$cellNames %in% keep]

# filter the ArchR object
filt_proj <- filt_proj[overlap, ]

# subset both the processed atac data and rna data to get only the overlapping cells
rna_overlap <- proc_rna[, colnames(proc_rna) %in% overlap] 

atac_overlap <- proc_atac[, colnames(proc_atac) %in% overlap]

Save the subset ArchR project and proceed with it for now. It might be interesting to revisit and see if I can cluster and identify the cell types de novo at a later stage!

Dimensionality reduction

• two other algorithms:
  • latent semantic indexing (LSI) in Signac
  • landmark diffusion maps (LDM) in SnapATAC
• ArchR: optimized iterative LSI method
  • exhibits less susceptibility to batch effects
  • focuses on most variable features
  1. create a LSI Reduction from a subset of the total cells
  2. linearly project the remaining cells into this subspace with LSI projection (based on SVD)

Because we can have maximally two accessible alleles per cell, the scATAC-seq data is sparse. Therefore, the majority of accessible regions are not transposed, meaning that most loci will have 0 accessible alleles. A zero could mean “non-accessible” or “not sampled”. For many analysis we can use a binarized matix. Imporantly, the 1s have information, BUT the 0s do not!

A PCA would result in high inter-cell similarity at all 0 positions. An alternative approach for dimensionality reduction is a layered dimensionality reduction. First, Latent Semantic Indexing (LSI) is used. LSI is an approach from language processing. Different samples are the “documents” and different regions/peaks are the “words”.

Iterative LSI

1. compute term frequency (depth normalization to 10,000 per single cell)

\(TF = \frac{C_{ij}}{F_{j}}\), the term frequency of \(C_{ij}\) which is tiles i in cell j and \(F_{j}\) being the total number of tiles with insertions in cell j.

2. Inverse document frequency

• weights features by how often they occur
• more weight to less frequent peaks

\(IDF_{i} =\log {\frac{N}{1 + DF_{i}}\) with N being the total number of cells in the dataset and \(DF_{i}\) being the total number of cells (documents) in which the tile i is present.

3. The term frequency TF is normalized by the inverse document frequncy IDF. You get a TF-IDF matrix (term frequency-inverse document frequency) which tells us how important a region/peak is to a sample. In other words you transform a binary matrix to a non-binary matrix.

\(TF-IDF = \log{1 + TF{ij}} * IDF_{i} * 10^{4}}\)

4. SVD identifies the most valuable information across samples. Then we can use these most valuable features to represent the data in a lower dimensional space
5. Clusters are identified with Seurat’s Shared Nearest Neighbor clustering
6. Sum accessibility across all single cells in each cluster -> log-normalize
7. Identify most variable features across the clusters
8. repeat with most variable peaks as features

With LSI we can reduce the dimensionality of the sparse insertion matrix to tens or hundreds. Then UMAP or t-SNE can be used to visualize the data

Unlike in scRNA-seq we cannot select the top highly variable features before dimensionality reduction (high noise, low reproducibility). Rather the iterative LSI approach first computes a LSI on the most accessible tiles (this will identify clusters corresponding to the major cell types). Then, ArchR computes the average accessibility across these clusters across all features. Next, the most variable peaks across these clusters are identified. The most highly accessible peaks are the features of a new round of LSI. We can set how many rounds of LSI we want to be peformed.

Using iterative LSI reduces batch effects. If you see some batch effects you could try to add more LSI iterations and start from a lower initial clustering resolution. Also, the number of variable features can be lowered. #

#filt_proj <- addIterativeLSI(ArchRProj = filt_proj, useMatrix = "TileMatrix", 
                           #  name = "atac_LSI_25000")
#filt_proj <- addIterativeLSI(ArchRProj = filt_proj, useMatrix = "TileMatrix", 
                          #   name = "atac_LSI_50000")
filt_proj <- addIterativeLSI(ArchRProj = filt_proj, useMatrix = "TileMatrix", 
                             name = "atac_LSI_100000")

Clustering

Calling clusters in this new space uses the Seurat’s graph clustering function as default clustering method. The Seurat method first computs KNN graph and then a modularity optimization technique to cluster the cells (iteratively group cells together with Louvian algorithm using 10 random starts). Another option is to use “Scran”. The default number of nearest neighbors used is 10. The minimum number of cells for a cluster to be called a cluster is set to 5 by default. The maximum number of clusters to be called is set to 25 by default.

proj <- addClusters(input = proj, reducedDims = "IterativeLSI")

Visualization

We will add the celltype annotation from the processed data to our ArchR project. Since we already have annotations, we do not need to perform clustering.

#filt_proj <- addTSNE(ArchRProj =filt_proj, reducedDims = "atac_LSI_100000")

filt_proj <- addUMAP(ArchRProj =filt_proj, reducedDims = "atac_LSI_100000")

filt_proj <- loadArchRProject(path = "02_filtered_annotated_ArchRobject/")


# add celltype annotations form the processed atac to the ArchR object

# get vector of celltype annotations
celltypes <- as.character(atac_overlap$celltype.mapped)

# name celltype vector with cell names
names(celltypes) <- colnames(atac_overlap)

#filt_proj <- addCellColData(filt_proj,
 #              data = celltypes, 
  #             name = "celltypes",
   #            cells = colnames(atac_overlap)
    #           )

# add custom colors
colPalette_celltypes = c('#532C8A',
 '#c19f70',
 '#f9decf',
 '#c9a997',
 '#B51D8D',
 '#3F84AA',
 '#9e6762',
 '#354E23',
 '#F397C0',
 '#ff891c',
 '#635547',
 '#C72228',
 '#f79083',
 '#EF4E22',
 '#989898',
 '#7F6874',
 '#8870ad',
 '#647a4f',
 '#EF5A9D',
 '#FBBE92',
 '#139992',
 '#cc7818',
 '#DFCDE4',
 '#8EC792',
 '#C594BF',
 '#C3C388',
 '#0F4A9C',
 '#FACB12',
 '#8DB5CE',
 '#1A1A1A',
 '#C9EBFB',
 '#DABE99',
 '#65A83E',
 '#005579',
 '#CDE088',
 '#f7f79e',
 '#F6BFCB')

celltypes <- (as.data.frame(getCellColData(filt_proj)) %>% group_by(celltypes) %>% 
 summarise(n = n()))$celltypes

col <- setNames(colPalette_celltypes, celltypes)

df <- as_data_frame(cbind(getCellColData(filt_proj), getEmbedding(filt_proj)) ) %>%
  rename(c(umap1 = atac_LSI_100000.UMAP_Dimension_1, umap2  = atac_LSI_100000.UMAP_Dimension_2))

variables <- c("celltypes", "Sample", "nFrags", "DoubletScore")


plot1 <- map(c("celltypes"), function(n){
  ggplot() +
  geom_point(aes(x = df %>% pull("umap1"), y = df %>% pull("umap2"), 
                 col = df %>% pull(n)),  size = .8) +
    guides(col=guide_legend(title=paste0(n))) +
    scale_color_manual(values = col)+
    xlab("umap1") +
    ylab("umpa2") +
    labs(title = paste0(n))
})

plots2 <- map(c("nFrags", "DoubletScore"), function(n){
  ggplot() +
  geom_point(aes(x = df %>% pull("umap1"), y = df %>% pull("umap2"), 
                 col = df %>% pull(n)), size = .8) +
    scale_color_viridis_c() +
    guides(fill=guide_legend(title=paste0(n))) +
    xlab("umap1") +
    ylab("umpa2") +
    labs(title = paste0(n))
})

plot3 <- map(c("Sample"), function(n){
  ggplot() +
  geom_point(aes(x = df %>% pull("umap1"), y = df %>% pull("umap2"), 
                 col = df %>% pull(n)),  size = .8) +
    guides(col=guide_legend(title=paste0(n))) +
    xlab("umap1") +
    ylab("umpa2") +
    labs(title = paste0(n))
})


plot1
## [[1]]

do.call(gridExtra::grid.arrange, c(plot3, plots2, ncol=3))#, nrow = 2))

#saveArchRProject(filt_proj, outputDirectory = "02_filtered_annotated_ArchRobject", load = FALSE)

Check gene scores for the different celltypes

filt_proj <- addImputeWeights(filt_proj, reducedDims = "atac_LSI_100000")
## Warning in sprintf("Completed Getting Magic Weights!",
## round(object.size(weightList)/10^9, : one argument not used by format 'Completed
## Getting Magic Weights!'

markerGenes  <- c(
    "Dnmt3b",  # Epigblas
    "Ptn", # Forebrain/Midbrain/Hindbrain
    "Pou5f1", # Primitive Streak
    "Myl7", # cardiomyocytes
    "Hoxa10" #Allantois

  )

#{r, fig.height=10, fig.width=10} p <- plotEmbedding( ArchRProj = filt_proj, colorBy = "GeneScoreMatrix", name = markerGenes, embedding = "UMAP", imputeWeights = getImputeWeights(proj) ) do.call(gridExtra::grid.arrange, c(p, ncol = 3))

Visualizing Genome Browser Tracks

Browse local chromatin accessibility at marker genes. Plot genome browser tracks per celltype

p <- plotBrowserTrack(
    ArchRProj = filt_proj, 
    groupBy = "celltypes", 
    geneSymbol = markerGenes, # the plot window is centered at the TSS
    upstream = 50000,
   
     downstream = 50000
)

stages_plots <- plotBrowserTrack(
    ArchRProj = filt_proj, 
    groupBy = "Sample", 
    geneSymbol = markerGenes, 
    upstream = 50000,
   
     downstream = 50000
)

markerGenes  <- c(
    "Dnmt3b",  # Epigblas
    "Ptn", # Forebrain/Midbrain/Hindbrain
    "Pou5f1", # Primitive Streak
    "Myl7", # cardiomyocytes
    "Hoxa10" #Allantois

  )

• Dnmt3b = Epiblast
• Ptn = Forebrain/Midbrain/Hindbrain
• Pouf5f1 = Primitive Streak
• Myl7 = cardiomyocytes
• Hoxa10 = Allantois

Since we group the cells accoring to celltype below and plot the average signal/accessibility. Since some signals are only expected in a subset of the embryonic stages, the signals might be diluted if we look at all stages at once and take the mean accessibility over different stages. It would probabl

```#{r,fig.height=15, fig.width=10} for (elem in p){ grid::grid.newpage() grid::grid.draw(elem) }

```#{r}
for (elem in stages_plots){
  grid::grid.newpage()
  grid::grid.draw(elem)
}

Different Time points

To differentiate between expression at different stages we will have to split the ArchR projects into

Combine with RNA-seq data

Since we are using multiome data, we do not need to integrate the ATAC and RNA data, but we have both informations from the same cell. We will use the processed scRNA-seq data which we filtered to contain only the cells which overlap with the cells in the ArchR object.

Look here for more details on multiome analysis in ArchR: https://greenleaflab.github.io/ArchR_2020/Ex-Analyze-Multiome.html

One interesting thing to notice is that the number of genes in the gene score matrix is different from the number of genes in the RNA dataset.

I get an error when trying to extract the tile matrix. Not sure what is happening there.

# extract gene score matrix
gene_scores <- getMatrixFromProject(filt_proj, useMatrix = "GeneScoreMatrix")
gene_scores_mtx <- assays(gene_scores)[[1]]

# set names of the gene score matrix
rownames(gene_scores_mtx) <- rowData(gene_scores)$name

#tile_matrix <- getMatrixFromProject(filt_proj, useMatrix = "TileMatrix", binarize = TRUE)

print(paste0("The number of genes in the ATAC-seq Gene Score Matrix is: ",
             nrow(gene_scores_mtx)))
## [1] "The number of genes in the ATAC-seq Gene Score Matrix is: 24333"

print(paste0("The number of genes in the RNA-seq Gene Expression Matrix is: ",
             dim(rna_overlap)[1]))
## [1] "The number of genes in the RNA-seq Gene Expression Matrix is: 32245"

print(paste0("The number of genes found in both the Gene Expression Matrix and the Gene Score Matrix:", 
      length(intersect(rownames(gene_scores_mtx), rownames(rna_overlap)))))
## [1] "The number of genes found in both the Gene Expression Matrix and the Gene Score Matrix:22273"

Filter genes which are found in gene score matrix and gene expression matrix. Reorder cell names to match between gene scores and gene expression. Here it would be nice to show VennDiagrams to quickly be able to understand what the differences in genes are between ArchR object and gene expression matrix. Is there a way to include also genes which do not overlap?

# Try to reorder the cellnames to the 

# this is the order of the cell names in the ArchR project gene score matrix
order <- rownames(getCellColData(filt_proj))

# extract the gene expression matrix 
rna_counts <- assays(rna_overlap)[[1]]

# reorder gene expression matrix cell names to match the cellnames of the ArchR object
rna_counts <- rna_counts[, order]

# I checked that they are the same
#colnames(rna_counts) == rownames(getCellColData(filt_proj))

# get the gene names of the gene score matrix
genes <- rownames(gene_scores_mtx)

# subset the gene expression matrix to obtain only the genes found in the gene 
# score matrix
rna_counts <- rna_counts[rownames(rna_counts) %in% genes,]

Next we have to create a GRanges object which will be used in the next step to create a RangedSummarizedExperiment. When using the gene annotation metadata stored in the rowData of the RNA-SCE for creating this RangedSummarizedExperiment with the ArchRobject, I got an error message. This might have been due to the fact that the gene annotations looked different between the RNA-SCE and the ArchR object. Therefore, I tried using the same gene annotation as ArchR uses. In order for this to work, we will have to exclude any genes which are not used by ArchR.

# get the Granges object from the Archr project
archr_granges <-getGenes(filt_proj)

# filter for genes which are included in the gene expression matrix
granges <- subset(archr_granges, symbol %in% rownames(rna_counts))

# the problem here was that the gene annotations of the expression matrix
did not match the gene annotations of the Archr object.

df <- rowData(rna_overlap) %>% as.data.frame() %>% filter(gene %in%genes)
df

granges <- df %>% GRanges()
granges

We will now add the gene expression matrix to the ArchR object. First, we have to convert the Single Cell Experiment to a Summarized Experiment.

# create summarized experiment of the scRNA-seq
seRNA <- SummarizedExperiment(assays = list(counts=rna_counts),
                          colData = colnames(rna_overlap),
                          rowRanges = granges)
seRNA
## class: RangedSummarizedExperiment 
## dim: 22273 28154 
## metadata(0):
## assays(1): counts
## rownames(22273): Xkr4 Rp1 ... Asmt Ccl21c
## rowData names(2): gene_id symbol
## colnames(28154): E7.5_rep2#TGTAATGTCATGCTAA-1
##   E7.5_rep2#GCGCTTAAGCATTTCT-1 ... E8.0_rep2#CGCCTGTGTACTGAAT-1
##   E8.0_rep2#TTGTGAGGTTTAAAGC-1
## colData names(1): X


# add RNA gene expression matrix to the ArchR object
filt_proj <- addGeneExpressionMatrix(input = filt_proj, seRNA = seRNA, force = TRUE)


# save the ArchR object with the added gene expression matrix
#saveArchRProject(filt_proj, outputDirectory = "03_added_gene_expression_ArchRobject", load = FALSE)

We will add a second LSI for the RNA data. The one based on the ATAC data is called “LSI_ATAC”. We will call the one for the RNA data “LSI_RNA”

filt_proj <- addIterativeLSI(
    ArchRProj = filt_proj, 
    clusterParams = list(
      resolution = 0.2, 
      sampleCells = 10000,
      n.start = 10
    ),
    saveIterations = FALSE,
    useMatrix = "GeneExpressionMatrix", 
    depthCol = "Gex_nUMI",
    varFeatures = 2500,
    firstSelection = "variable",
    binarize = FALSE,
    name = "LSI_RNA"
)

filt_proj <- addIterativeLSI(
  ArchRProj = filt_proj, 
  saveIterations = FALSE,
  clusterParams  = list(
    resolution = 0.2,
    sampleCells = 10000,
    n.start = 10
  ),
  useMatrix = "GeneExpressionMatrix",
  depthCol = "Gex_nUMI",
  varFeatures = 100000, 
  firstSelection = "Var", 
  binarize = FALSE,
  name = "LSI_RNA"
)

Now we will combine the two embeddings.

#Combined Dims
proj <- addCombinedDims(proj, reducedDims = c("LSI_ATAC", "LSI_RNA"), name =  "LSI_Combined")

#UMAPs
proj <- addUMAP(proj, reducedDims = "LSI_ATAC", name = "UMAP_ATAC", minDist = 0.8, force = TRUE)

proj <- addUMAP(proj, reducedDims = "LSI_RNA", name = "UMAP_RNA", minDist = 0.8, force = TRUE)

proj <- addUMAP(proj, reducedDims = "LSI_Combined", name = "UMAP_Combined", minDist = 0.8, force = TRUE)